package scu.maqiang.heat;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.mesh.T3Type;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class StaticRadiation {

	public static void main(String[] args) {
		Mesh2T3 mesh = new Mesh2T3().square2D(20, 20, T3Type.Left);
		FES2T31 fs = new FES2T31(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		fs.assembleStiff(new double[]{4.0}, BVPType.COMMON, A);
		double[] massCoef = new double[fs.GetNdof()];
		Arrays.fill(massCoef, 1.0);
		double[][] mCoef = new double[][] {massCoef};
		SRMatrix M = new SRMatrix(fs.GetNdof());
		
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleSource(new double[]{1.059300213000000e+03}, BVPType.COMMON, RHS);
		NewIterSSolver solver = new NewIterSSolver();
		solver.setEps(1.0e-6);
		double[] x0 = new double[fs.GetNdof()];
		double[] x = new double[fs.GetNdof()];
		fs.applyBC_MBN(RHS, Direct.All, 1000, 1);
		while(true) {
			fs.assembleMass(mCoef, BVPType.CONSIST_MASS, M);
			SRMatrix LHS = A.add(1.0, M);
			fs.applyBC_MBN(LHS, Direct.All, 1);
			solver.setMatrix(LHS);
			solver.PCGSSOR(RHS, x, 1.5, 1);
			if(MVO.distance(x0, x) < 1.0e-10) {
				break;
			}
			Arrays.setAll(x0, i -> x[i]);
			Arrays.setAll(massCoef, i -> 2.0 + 5.670373e-08 * x[i] * x[i] * x[i]);
			LHS.clear();
			M.clear();
		}
		mesh.toTecplot("Radiation.dat", x);
	}

}
